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An  Algorithm  For  Computing  Matrix  Square  Roots  With  Application  To  Riccati  Equation  Implementation 


D.  W.  Repperger 


Aerospace  Medical  Research  Laboratory 
Wright-Patterson  Air  Force  Base,  Ohio  45433 


Abstract 


An  iterative  algorithm  is  presented  for  obtain- 
ing a positive  definite  symmetric  square  root  of  a 
positive  definite  symmetric  matrix.  This  algorithm 
has  application  in  the  implementation  of  Riccati 
type  equations.  The  approach  presented  here  has  the 
advantage  that  apriori  upper  and  lower  bounds  to 
the  converged  answer  can  be  obtained  sequentially 
at  any  point  in  the  iterative  process.  These 
apriori  bounds  can  also  be  obtained  for  the  Riccati 
equation  using  a discrete  implementation  procedure. 
Theorems  on  convergence  are  proved  and  examples 
are  worked. 


I.  Introduction 


The  study  of  optimal  control  and  estimation  the- 
ory has  been  influenced  in  recent  years  by  factori- 
zation methods  in  the  analytic  expression  of  the 
filtering  equations  [1,2,3].  Computationally  the 
solution  of  filtering  and  smoothing  equations  appears 
numerically  better  behaved  [4]  when  the  equations 
are  of  the  square  root  type  in  lieu  of  the  standard 
forms.  In  the  study  of  Riccati  type  equations  [5], 
an  attempt  has  been  made  to  apply  the  factorization 
methods  directly  to  the  steady  state  Riccati  equa- 
tion but  an  algorithm  to  produce  these  results  is 
required. 

In  the  study  of  Riccati  type  equations,  many  meth- 
ods exist  to  determine  solutions  such  as  iterative 
procedures  [6,7],  the  partitioned  algorithm  approach 
[8,9,10],  the  Matrix  Sign  Function  method  [11,12], 
and  other  methods  applicable  in  numerical  integra- 
tion, e.g.  [15].  The  Matrix  Sign  Function  approach 
is  one  method  which  has  the  advantage  of  obtaining 
additional  non-positive  definite  solutions  through 
the  use  of  a Symplectic  matrix  composed  of  the  ma- 
trices in  the  Riccati  equation.  This  aspect  has 
been  discussed  by  Bucy  [13]  and  Potter  [14]. 

This  paper  will  develop  an  algorithm  for  deter- 
minlng  the  square  root  of  a positive  definite 
*The  research  reported  in  this  paper  was  sponsored 
by  Aerospace  Medical  Research  Laboratory,  Aero- 
space Medical  Division,  Air  Force  Systems  Command, 
Wright-Patterson  Air  Force  Base,  Ohio  45433.  Fur- 
ther reproduction  is  authorized  to  satisfy  needs 
of  the  U.S.  Government. 


symmetric  matrix.  The  convergence  of  the  algorithm 
will  be  proved  with  upper  and  lower  apriori  bounds 
on  the  converged  value  determined  at  any  point  along 
the  sequential  process.  This  algorithm  can  then  be 
applied  to  the  implementation  of  Riccati  type  equa- 
tions. The  Riccati  type  equations  that  are  applicable 
for  this  approach  include,  but  are  not  limited  to, 
the  discrete  version  of  the  partitioned  Riccati 
solutions  [8,10],  and  the  various  covariance  express- 
ions which  occur  in  square  root  filtering  methods 
already  discussed  [3,15].  It  is  noted  that  the  square 
root  algorithm  presented  here  will  be  compared 
(through  examples)  to  the  Matrix  Sign  approach  [11] 
and  not  to  the  triangular  factorization  methods.  In 
a hard  to  find  reference  [16]  another  (but  different) 
algorithm  was  developed  which  has  been  discussed  by 
Bellman  [17].  The  algorithm  in  [16],  however,  has 
limitations  on  the  norm  of  the  matrix  considered  and 
appears  to  be  related  to  a version  of  the  spectral 
factorization  problem.  In  Astrom's  book  [18],  similar 
algorithms  result  from  studies  of  related  spectral 
factorization  problems. 

In  order  to  better  understand  why  the  factoriza- 
tion methods  yield  algorithms  which  have  numerical 
advantages  over  other  methods,  this  paper  is  divided 
in  four  parts.  First  the  motivation  for  using  this 
approach  is  demonstrated  by  working  a scalar  example 
using  an  algorithm  from  Number  Theory  called  Euclid's 
algorithm  to  characterize  an  irrational  number  in 
terms  of  continued  fractions  . Such  methods  are  used 
in  Number  Theory  [19,20],  and  for  scalar  irrationals 
(such  as  a square  root),  they  give  rise  to  definitions 
such  as  "most  efficient  expansion"  and  "best  possible 
approximation".  By  defining  a structure  of  the  contin- 
ued fraction  expansion  subject  to  certain  constraints, 
the  definition  of  "most  efficient  expansion"  can  be 
given  in  an  explicit  manner.  Using  a sequential  proce- 
dure defined  as  Euclid's  algorithm  [19,20],  the  contin- 
ued fraction  expansion  is  obtained  for  the  irrational 
number  it  . Euclid's  method  ( or  a "most  efficient  ex- 
pansion") is  then  constructed  for  a scalar  square  root. 
Apriori  upper  and  lower  bounds  are  obtained  for  both 
expansions  of  these  two  scalar  irrational  numbers. 

Part  II  of  this  paper  introduces  the  "square  root 
algorithm"  which  differs  from  Euclid's  method.  The 
square  root  algorithm  is  applied  first  to  the  scalar 
square  root  and  the  resulting  continued  fraction  ex- 
pansion is  compared  to  Euclid's  method.  The  apriori 
upper  and  lower  bounds  are  also  calculated.  The  third 
part  of  this  paper  extends  the  scalar  square  root  al- 
gorithm to  the  matrix  case.  Theorems  on  convergence 
are  proved  and  examples  are  worked.  The  matrix  square 
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root  algorithm  considered  here  is  compared  (through 
an  example)  to  the  Matrix  Sign  Approach  and  the  con- 
vergence of  the  algorithm  is  studied  numerically. 

Finally,  the  fourth  part  of  this  paper  considers 
the  application  of  this  approach  to  the  implementa- 
tion of  Riccati  type  equations.  A manner  of  obtain- 
ing apriori  upper  and  lower  bounds  is  demonstrated 
for  Riccati  type  equations.  Some  theorems  are  given 
and  an  algorithm  is  demonstrated  which  extends  the 
results  obtained  here  for  the  matrix  case  into  an 
algorithm  for  computing  the  Riccati  equation.  The 
apriori  upper  and  lower  bounds  for  the  Riccati  equ- 
tion  are  obtained  sequentially  as  in  the  matrix  case, 
First,  part  I of  this  paper  will  introduce  Euclid’s 
algorithm  and  provide  a motivation  for  using  fac- 
torization algorithms. 


II.  Part-I  Euclid’s  Algorithm-A  Motivation  For 
Using  Square  Root  Methods 


It  is  well  known  that  numerical  problems  occur  in 
computing  Riccati  solutions  or  other  matrix  equations 
if  the  matrices  are  ill  conditioned  [21].  One  measure 
of  ill  conditioning  occurs  if  the  eigenvalues  of  the 
free  system  are  seperated  by  more  than  one  order  of 
magnitude.  One  would  expect,  therefore,  that  factor- 
ization methods  should  prove  to  have  numerical  advan- 
tages because  the  square  root  matrix  may  have  eigen- 
values with  less  dispersion  then  the  original  free 
system  matrix.  In  an  effort  to  study  this  effect  num- 
erically, Euclid's  algorithm  is  introduced  and  a sca- 
lar example  is  worked  from  Number  Theory  on  an  irra- 
tional number.  The  purpose  of  this  example  is  to  ill- 
ustrate to  the  reader  the  relationship  between  the 
factorization  methods  to  continued  fraction  approaches 
and  also  to  define  the  numerical  advantages  (in  terms 
of  a definition  called  "most  efficient")  of  the  fact- 
orization methods  considered  here.  First  some  defini- 
tions are  necessary. 

Let  denote  a scalar  irrational  dumber.  The  pro- 
blem of  interest  is  to  determine  a continued  fraction 
expansion  of  in  terms  of  rational  numbers.  First 
the  definition  of  a continued  fraction  is  specified 
as  follows  [19,20]: 

Def ine: 

w “ <?o>  al'  •••  * aN>  ’’>>  d) 

which  is  notation  for: 

1 

W*an-f  a,  + 1 

a2  + 1 " ~ 

a 3 + 

* * + _1 

*N  + • 


•(2) 

where  each  a^  is  subject  to  the  constraint  of  being 
positive  and  an  Integer.  Then  Euclid's  algorithm  [20] 
is  defined  as  follows: 

Denote:  a±  =■  [5  i]  (3) 

where  a^  is  defined  as  the  nearest  integer  smaller 
than  the  irrational  number  Euclid's  algorithm  for 
computing  the  a^  from  the  ^ proceeds  as  follows: 


1 

£1  " So  - ao 

and  let  = [€  ^].  Inductively  it  follows  that: 


and  a.  = [5^  ] and  this  completes  the  algorithm  for 
all  a1.  It  is  now  appropriate  to  give  a definition  of 
"most1ef ficient  expansion"  as  it  relates  in  the  con- 
text of  the  continued  fractions  in  equations  (1,2). 
Definition:  "Most  Efficient  Expansion" 

The  most  efficient  [19,20]  expansion  of  an  irra- 
tional scalar  number  by  a continued  fraction  rep- 
resentation is  the  one  in  which  E0  can  be  accurately 
expressed  to  as  many  decimals  as  possible  by  the 
fewest  number  of  integer  terms  in  the  expansion.  It 
is  noted  that  W is  subject  to  the  constraint  of  equa- 
tion (2)  with  each  a.  positive  integers.  It  has 
been  shown  [19,20]  that  Euclid's  method  satisfies 
this  property  for  scalar  irrationals.  An  example 
will  now  be  worked  with  the  irrational  number  * 
to  illustrate  Euclid's  algorithm. 

Example  (1) : 

Let  be  the  irrational  number  it,  i.e.  - 3.14159+ 
Hence  a0=[E  o]  = 3-  T°  calculate  a^,  compute 


or  1 1 

6 . = C - 3 = .14159  = 7.06+ 

1 0 

Now  compute  « [5  ^]  = 7.  In  this  manner  we  calcu- 
late (using  double  precision  (29  digits)  on  a 
CDC  computer)  the  following  numerical  results: 


TT 


e 5 
« 6 


5 8 

E 9 
5 10 

e 11 

£ 12 


3.141592653589793238462643383279 

7.0625133059310457697930051531 

15.99659440668569411310599874 

1.0034172310133977856369454934 

292 . 63459101223866070785852178 

1.5758180949841629954527421461 

1. 7366595609113341887024369765 

1.3574791573503535151332712575 

2.7973658867611542733590927681 

1. 2541293985649846777332901666 

3.9350032135077287460064498245 

1.0695150407541735890218572566 

14.385376015764735321479897720 


This  gives  rise  to  aQ=3,  aj  = 7,  32=13,  aj=l,  a^=292, 
85=1 , a^=l,  a^— 1,  ag— 2,  aq*l,  a,Q  = 3,  a^^=l,  and 
aj2=l^  which  results  in  tne  following  partial 
fraction  expansion: 


n 


1 


3+  7+1 

13+_1 

1 +_1 

292+  1 

1+1 

1+1 

1+1 

2+1 

1+1 

3+  1 

1+  1 

14+ 


aQ  = [to],  now  proceeding  as  an  algorithm,  let, 


(A) 


i 
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In  order  to  study  convergence  and  to  develop  apriori 
bounds,  it  is  of  interest  to  study  the  partial  suras 
of  equation  (4).  Let  r0=a0,  r^=  a + 1/a^,  and  in 
general: 

CN  = <?°>  al*  •••  > ^ {5) 

It  has  been  shown  explicitly  [19,20]  that  for  the 
scalar  case,  the  following  upper  and  lower  bounds 
exist  for  tt  : 

ro<r2<r4<  •••<  rN  = " <""<r5<r3<rl 

lira  N-»-ao  (6) 

which  gives  rise  to  apriori  bounds  of  the  form: 

/ 333  /..../* <T...  355  / 22 

^ 106  ^ ^ ^ ^ 113  ^ 7 


1 . 3660254037844363479757026605 
2.7320508075688992137673133780 
1.3660254037843965712518405766 
2.7320508075691909097423020388 
1.3660254037838529560257256617 
2.7320508075732481355762416877 
1.3660254037762821195840238963 
2.7320508076297576012776604042 
1.3660254036708340246350351294 
2. 7320508084168328955065207969 
1.3660254022021315334827392008 


which  results  in  a =1 ,a, =1 ,a2=2 ,aj=l ,a^=2, a ^=1, 
a^=2)a2=l,ag=21aj=l,  and  a^g=2.  In  this  manner  the 
"most  efficient"  method  to  calculate  the  square 
root  of  3 is  as  follows  using  the  continued  fraction 
approach : 


The  scalar  version  of  this  proof  was  presented  for 
the  square  root  of  a scalar  number  in  [22]  , the 
matrix  case  will  be  presented  in  the  sequel.  There- 
fore by  calculating  rQ  and  r ^ , an  upper  and  lower 
bound  on  it  can  be  determined  apriori.  By  then  calcu- 
lating r^  and  r^,  an  even  more  accurate  bound  on  it 
can  then  be  determined.  This  procedure  can  be  con- 
tinued indefinitely  with  convergence  to  values  as 
close  as  desired.  To  illustrate  to  the  reader  how 
these  apriori  bounds  can  be  obtained  for  the  irra- 
tional number  x , the  terms  rQ,  ••••,  rg  are  calcu- 
lated for  the  fraction  in  equation  (4) . 


3 =3 
22/7  = 3 
333/106  * 3 
355/113  = 3 
103993/33102  = 
104348/33215  = 
208341/66317  = 
312689/99532  = 
833719/265381  ■ 


.142857142857+ 

,14150943396226415094+ 

.14159292035+ 

3.141592653011+ 

3.14159265392142+ 

3.141592653467+ 

3.1415926536189+ 

■ 3.14159265358107777+ 


1+  1 

2+  1 


The  periodic  nature  of  the  a^  terms  in  equation  (8) 
is  of  interest  because  it  is  known  that  in  Riccati 
equation  expansions  of  the  partitioned  algorithm 
approach  [8],  this  same  periodicity  occurs.  The 
periodicity  obtained  here  numerically  occurs  due  to: 

5 j = € j+2  for  a11  J 

with  a^  = 1 (i  odd) 

ai+1  = 2 (i  odd) 

and  this  is  a result  due  to  the  following  identity: 


It  is  easily  seen  from  these  numerical  values  that 
the  r.  satisfy  the  relationships  specified  by  equa- 
tion *6) . Also  from  the  example  it  is  observed  that 
the  first  four  terms  of  the  expression  r ^ give 
accuracy  beyond  7 decimal  places.  One  should  now  be 
motivated  to  apply  Euclid’s  algorithm  to  a quadratic 
equation  to  observe  the  results. 

Ill  Euclid’s  Algorithm  Applied  To  A Scalar  Square 
Root 

Taking  as  an  example  the  square  root  of  3,  the 
most  efficient  method  will  be  studied  numerically 
using  example  (2) : 

Example  (2) : 

Expanding  the  square  root  of  3 to  30  decimal  places 
yields : 


The  results  obtained  here  so  far  indicate  that  using 
the  "most  efficient"  approach,  the  square  root  meth- 
ods have  a periodic  property.  If  this  method  were  to 
be  applied  to  the  matrix  case,  the  results  do  not  ex- 
tend readily  due  to  the  fact  that  the  definition 
"most  efficient"  does  not  have  the  same  meaning  in 
the  matrix  case.  In  order  to  have  an  algoritlim  that 
does  extend  to  the  matrix  case,  the  square  root  al- 
gorithm will  next  be  presented.  The  square  root  al- 
gorithm gives  identical  results  to  the  "most  effi- 
cient" algorithm  for  the  scalar,  case.  The  square 
root  algorithm,  however,  can  be  extended  to  obtain 
matrix  square  roots,  to  implement  the  Riccati  equa- 
tion, and  , in  addition,  to  determine  upper  and  low- 
er apriori  bounds. 

IV  - Part  II-  The  Square  Root  Algorithm  (Scalar  Case) 


* /5”  = 1 . 732050807568877293527446341505.  Now  Consider  example  (2)  to  determine  the  square  root 

the  calculation  of  the  a^  and  will  proceed  as  for  of  3 via  the  following  sequence  of  steps:  — | 


1.732+  ] - 1 


i = VT  - i 


Rewrite  this  as 


Let  W = ^7 

W = 1 + ( \fJ  - 1)  (9) 


Hence  a.  = C “ 1,A  summary  of  the  first  11  terms  or  W » 1 + 1 

yields:  1 

jf  ( \|T  + l) 


l + l 

| + y(i+\TT  - i) 


.,11  COBB 
0 of  SPEtUL 


ft  \0 


amo  | 


W “ 1 + 1+|(<3  -1) 


1 + 1+lr 

2 1 

I ( 5 +i) 


W = 1 + 1 + 1 

2 + 1 

1 + 1 

2 + 1 

1 +1 
2+ 


Since  3 = 1.7320508+  , then  the  continued  fraction 

of  interest  is: 

1 

1.732+  = 1 + 1+  1 

2+  1 

1+  1 

2 + 1 

1 + 1 

2+ 


This  square  root  approach  yields  numerical  answers 
for  this  example  which  are  identical  to  the "most  eff- 
icient1' method.  For  completeness  the  partial  sums  r0, 
r^,  ate  computed  and  displayed  here: 


to  = 1 


= 

2/1  = 

2.0 

= 

5/3  = 

1.667 

= 

7/4  - 

1.750 

= 

19/11  = 

1.7272 

= 

26/15  = 

1.7333 

= 

71/41  = 

1.7317 

= 

97/56  = 

1.73214 

It  is  easily 


= 265/153  = 1.732026 
= 362/209  = 1.732057 
= 989/571  = 1.732049 
demonstrated  that: 


where  the  inequality  representation  for  two  matrices 
Ai>B  implies  that  A-B  is  positive  definite  and  I is 
the  nxn  identity  matrix.  If  equation  (15)  is  not  sat- 
isfied for  a specified  SQ  and  S^,  then  we  require: 

SQ  - e I>0,  and  Sx  - e 2I  >0  (16) 

for  some  e >0.  The  ensuing  derivation  follows  exac- 
tly with  (16)  replaced  by  (15)  , but  for  notation 
simplicity,  the  inequality  of  equation  (15)  will  be 
used.  To  derive  the  square  root  algorithm,  proceed 
as  in  equations  (9-12)  for  the  scalar  case. 

Denote:  , 


This  can  be  written  as: 


Theorem  1 will  now  be  helpful  in  the  derivation: 
Theorem  1: 

s*/2  - i = (Si-iKsJ^+i)-1  - (s|/2+i)_1(s1-i) 

sJ/2  + I =(S1-I)(S2/2-I)“1  - (S2/2-I)_1(S1-I)  (19b) 

Proof : Consider  the  identities: 

(S2/2  -I)(S2/2+I)  - Sj-I  (20a) 

(Sy2+  I)(Sy2-I)  - S.-I  (20b) 


(S2/2+  I)(S2/2-I)  ■=  Sx-I 


ro<r2<r4<r6<r8<r10  <•  • -MT<-  • -<r9 <t7<r5<r3 <1^ 

(13) 

and  approximations  to  the  square  root  of  3 can  be  det- 
ermined by  calculating  only  rQ  and  r^;  or  more  accu- 
rately be  calculating  and  r^,  etc.  This  approach 
will  now  be  extended  to  the  matrix  case. 

V-Part  Ill-Extension  of  The  Algorithm  For  Determining 
Matrix  Square  Roots 


Equation  (19a)  follows  by  pre  and  post  multiplying 
equations  (20a-b)  by  (S^/2  4-1)  1 which  is  positive 
definite.  Equation  (19b  ) follows  by  pre  and  post 
multiplying  equations  (20a-b)  by  (S^**  -i)“l  which 
is  known  to  be  positive  definite  by  the  inequality  (15). 

Q.E.D. 

The  derivation  of  the  algorithm  will  now  continue  by 
rewriting  equation  (18)  as  follows: 

S = I + (S2/2-I)(sy2+I)(S2/2+I)'1  (21) 

O 1 -L  1 

From  equation  (20a)  this  implies  (21)  can  be  written: 
s0  - I + (S1-I)  [S2/2  + I]  _1  (22) 


Since  for  two  matrices  A,B,  (AB) _-*-=B_^A  , it  is 

desired  to  represent  equation  (22)  as: 

sQ  - 1 + 1 [ (S1/2  +i)(s1-i)-1  r1  < 

This  can  be  written  as: 


In  order  to  develop  an  algorithm  for  matrix  square 
roots,  it  is  necessary  to  accurately  define  the  ma- 
trices of  interest.  Let  be  an  (nxn)  positive  def- 
inite symmetric  matrix.  The  (nxn)  matrix  SG  which 
is  the  square  root  of  is  required  to  be  positive 
definite  symmetric  and  satisfy: 


where  T denotes  matrix  transpose.  S is  defined  as  the 
the  square  root  matrix  of  S^.  For  simplicity  in  the 
ensuing  derivation,  it  will  be  assumed  that: 


SQ  - 1 > 0 and  Sj  - I > 0 


SQ  - I + I i ( S1-I)_±  + S1  (S1-I)  * ]“A  (24) 

Now  the  procedure  used  in  equations  (17-18)  will  be 
reapplied  to  equation  (24)  resulting  in: 

So  = I + X [ (Si-I)'1  + [I+sJ/2-I](S1-I)-1r1  (25) 

It  is  worthwhile  to  rewrite  this  as: 

Sc  - I + I (2(S1-I)'1  + [ S;/2  -I](S1-I)'1  r1 
Using  the  results  of  equation  (19b)  yields: 

S0  - I + I [ 2 (Sj-I)-1  +[S;/2  + I]-1  l'1  (26) 
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r 


k 


Now  expanding  the  inner  term  yields: 
So=I+I[  2(S1-D'1  +I[  I+(s}/2  -*-1-1) J-1  ]' 


(27) 


This  can  be  rewritten  as: 


X + I[  2(S1-I)"1  + I[2I+(S2/2-I)]_1  )-1 


(28) 


.1/2 


The  observant  reader  may  now  compare  the  term  (S^  -I) 

in  equation  (18)  with  the  same  term  in  the  inner  most 
brackets  of  equation  (28).  The  next  step  in  the  der- 
ivation is  to  repeat  the  procedure  from  equations 
(18-28)  resulting  in  the  next  expression  as  follows: 


S0=I+I[2(S  -I)'1+I[2I+(S,/2-I)(s}/2+I)(sy2+I)-1r1r1 

1 (29) 

etc.  Since  the  derivation  is  periodic  from  this  point 
on,  the  results  can  be  summarized  as  follows: 


Let 


SA  - 21 


2(S1-X) 


-1 


(30a) 

(30b) 


Then : 


S = I+I[SR+I[S.+I[S n+I[SA+I[S  *•*]  ]"A]  A] 


• • i _1i — 1 1 — 1 1 ~ i 


This  could  also  be  written  as  follows: 
I 


(31) 


S = I + S„  + I 
o B 


SA  + 1 


SB+  1 


SA  + 1 


SB  + 


The  expressions  (31-32)  obtained  are  desired  for 
several  reasons.  First,  the  matrix  case  reduces  to 


(32) 


the  "most  efficient"  approach  as  and  S become 
scalars.  Secondly,  the  matrices  S^  and  SgB  are 
positive  definite  and  it  is  desireable  to  have  this 


property  to  obtain  apriori  bounds.  Thirdly,  it  ap- 
pears that  this  method  has  application  in  the  im- 
plementation of  the  Riccati  equation  as  observed  in 
references  [3,8,15].  The  theorems  on  convergence 
will  now  be  given  with  their  proofs  in  Appendix  (A). 
Some  examples  will  then  be  worked  to  illustrate 
this  approach. 


Vl-Convergence  Proofs  of  The  Square  Root  Algorithm 


In  order  to  establish  a theorem  on  the  convergence 
of  the  nested  matrix  sequence  (31),  a sequence  of 
partial  sums  of  the  expression  (31)  will  be  derived. 
Denote  R^,  i*0,l,#,#  as  the  matrix  partial  sums  of 
(31)  which  reduce  to  the  scalar  partial  sums  r^ 
defined  in  equation  (5). 

Let : 

Ro  - 1 , (33a) 

RX  = I + I[  2(SrI)  i]  1 (33b) 


= I + 1[2(SX-1)-1  +H2I]-1]-1  (33c) 


R3  - I+I(2(S1-I)_1+I[2I+(2(S1-I)"1]-1]"1]-1 


I + I(2(S1-l)'1]~1 


(34a) 

(34b) 


Bj,  - I +[2(S1-l)'1+I[I+RN_2f1]'i(34c) 


Theorem  2 contains  the  properities  of  convergence  of 
this  algor ithm. 


Theorem  2: 


.1/2 


(35a) 


lira  R„ 

N— ► 0D  ‘ 1 

with  the  following  apriori  bounds: 

1/2 

Ro<'R2<R4<R6<'"<RN  * S1  ^ ••■<R7^R5<R3^R1 


lira  N- 


•00 


(35b) 


The  proof  of  theorem  (2)  is  given  in  appendix  (A). 
The  scalar  version  of  this  proof  was  given  in  [22] 
and  similar  proofs  of  the  scalar  example  are  given 
in  the  references  on  Number  Theory  [19,20].  Some 
examples  will  now  be  worked  to  illustrate  the  con- 
vergence properities  of  the  algorithm  and  also  how 
the  apriori  bounds  may  be  used. 


VII  Some  Numerical  Examples  Using  The  Matrix 
Square  Root  Algorithm: 


As  the  simplist  example  to  observe  the  proce- 
dure presented  here,  consider  the  square  root  of 
the  following  matrix: 


Example  3 


4.0  , 0.0 


0.0  , 25.0  ! 


Using  the  algorithm  results  in  R^q  and  R^  as 
follows : 


Ricr 


fl. 999797,  0.0 


L_° 


.0 


, 4.746472 


and  Rxx= 


2.000068,  0.0  | 
0.0  , 5. 176476J 


Please  recall  the  inequality: 


1/2  . 
■s.  <si 


is  true  for  all  N even  integers  and  all  M odd  inte- 
gers. Also,  the  apriori  bounds  Rxq  and  R indicate 
the  accuracy  of  this  procedure. 


Example  4 


4.0 

6.0 


6.0 

25.0 


06) 


(33d)  using  the  Matrix  Square  Root  Algorithm  results  in 


R and  R^^  as  follows: 


or  as  an  algorithm: 


938 


r 1.788847 

, .894398 

Tl. 788860 

1 , .8944^7 

R30=  | 

! .894398 

, 4.919241 

■ R31=  1 

1 .894447 

, 4.919423 

1 

— 

1 

— 

It  is  interesting  to  compute  the  products 

R30RqnT 

and  ^31^31  • 

T 

H3. 999921371  , 5.999700491~| 

R R_  « : 

30  30 

5.999700491  , 24.9988798 

T 

~4. 0000555, 

6.0002036 

R R - , 

31  31 

1 6.0002036, 

25.00075808 

Example  6 


S 


1 


4.0 

10.0 


10.0 

25.0 


(37) 


The  algorithm  yields: 


.772480  , 

1.844994 

R ■ 

”.771536  , 

1.845506 

1.844994, 

4.  646968_ 

*31 

f 

1.845506, 

4.647099 

The  interesting  result  occurs  by  considering  the 
products  RjqRjq1  and  R^R^. 


which  can  be  compared  to  of  equation  (36).  The 
next  example  was  taken  from  reference  [11]  to  serve 
as  a comparison  to  the  approach  presented  here  and 
also  to  look  at  higher  order  systems. 

Example  5 _ 

1.2  x 105  , 230.  , 10. 

S1  = 230.  , 1000.  , 1. 

10.  , 1.  , 0.5 

Using  the  Matrix  Sign  Approach  [11] , Denman  found 
the  computed  square  root  as: 

|346.40961  , 0.608420  , .0287555~ 

I 

S = I 0.608420  , 31.61690  , .0303966 

o , 

i 0.0287555  , 0.0303966,  0.705876 


R30R30 


4.00072821  , 9.998849043 

9.998849043  , 24.998314450 


R 


T 

31R31 


4.001202784 

10.000123404 


10.000123404 

25.001421516 

—J 


which  can  be  compared  to  S.  of  equation  (37);  this 
procedure  appears  to  work  if  S is  only  positive 
semi-definite.  Theorem  (3)  illustrates  that  a best 
estimate  of  a matrix  square  root  can  be  obtained 
by  taking  the  mean  value  of  the  best  two  apriori 
bounds . 

Theorem  3 : 

(3a)  The  "best  estimate"  of  a matrix  square  root 
can  be  obtained  as  follows: 


Using  the  Matrix  Square  Root  Algorithm  presented 
here,  the  results  obtained  were: 


SQ  = max  (1/2)  [ Rj,  + R^]  (38a) 


328.739943 
R632=  .574269 

I .027283 
'*364. 918997 

R633*  j • 644195 
. .030298 


.574269  , .027283 
31.616842,  ,030394| 

I 

.030394  , .705868J 

.644195,  . 030298  I 

l 

31.616978,  .030400 | 
.030400  , .7058681 


Obviously  for  higher  order  systems,  some  modes  con- 
verge much  faster  than  other  modes.  In  other  words 
the  ill  conditioning  of  will  effect  this  algor- 
ithm (by  slowing  down  its  convergence)  in  a manner 
dependent  on  the  degree  of  ill  conditioning  of  S3- 
Also  the  effects  of  numerical  bias  occur  in  this 
algorithm  as  can  be  seen  because  the  element  in  the 
third  row,  third  column  oA  Rg32  ant*  i,,  has  appar- 
ently converged  independent  of  the  remaining  elements. 
Continuing  the  calculation  >of  for  N>  2000  gave  no 
additional  change  in  this  term.  This  element  may  be 
numerically  biased  (through  roundoff  errors  or  trun- 
cation errors  involved  in  the  matrix  inversion  sub- 
routine used  here).  An  excellent  discussion  of  this 
problem  can  be  found  in  Bierman  [4].  One  more  ex- 
ample was  worked  to  numerically  study  S in  the  event 
it  was  only  positive  semi-definite. 


(3b)  The  error  in  the  estimate  specified  by  equa- 
tion (38a)  can  be  no  worse  than  the  following 
bound : 

II  So  -%W  * IIRIH-1-Rh|I  (38b) 

Hence  even  if  SQ  is  not  known.  Rjj,  ^ and  Rj^  provide 
an  estimate  of  the  accuracy  of  this  procedure.  One 
may  now  take  appropriate  matrix  norms  and  examine 
examples  3,4,  and  5.  Appendix  B demonstrates  a proof 
of  this  theorem  with  an  appropriate  definition  of 
"best  estimate".  This  result  can  be  seen  numerically 
in  examples  (3-5).  The  estimate  specified  by  equa- 
tion (38a)  is  as  accurate  a method  as  possible  to 
guess  a matrix  square  root.  The  methods  obtained 
here  will  now  be  applied  to  the  discrete  version  of 
the  Riccati  equation. 

VIII  Part  IV-Application  of  This  Approach  To 

Discrete  Implementation  of  The  Riccati  Equation 

In  the  implementation  of  the  continuous  Riccati 
equation  at  discrete  times  0=tQ<t^<t.,  with 

following  recursion  relationships 
[ W + P"1  J"1  AT  (39) 


A- 

can 


*kn  - v c 

be  obtained: 


the 


k+ 1 


P + A 


939 


where  PQ,  A,  and  W are  quantities  obtained  over  the 
time  intervals  of  a Riccati  equation  with  zero  in- 
itial conditions.  Equation  (39)  appears  in  one  form 
or  another  in  the  partitioned  algorithm  approach  of 
Lainiotis  [8],  the  square  root  filtering  algorithms 
of  Morf  and  Kailath  [3],  and  in  Potter  and  Womble 
[15].  To  apply  this  approach  similar  to  the  manner 
in  which  the  matrix  algorithm  (34a-c)  was  obtained, 
note  the  following  sequences: 


The  results  of  this  last  section  seem  most  interest- 
ing in  the  study  of  the  accuracy  of  methods  to  im- 
plement the  Riccati  equation.  It  is  the  contention 
here  that  since  these  methods  reduce  to  the  "most 
efficient"  method  in  the  scalar  case,  they  should  be 
considered  first  in  the  calculation  of  Riccati  type 
equations . 

IX  Summary  and  Conclusions 


- P 

o 

- Po  + A[W]_1AT 

= P +A[W+rPj-1]-1AT 
o ° 

- P0+A[W+[P0+A[W]'1AT]‘1]'1AT 
Algorithm  becomes: 

= P (40a) 

o 

- P + A[W]-1AT  (40b) 

O 


= Po  + A[W+(Rn_2)~1]_1AT  (40c) 

the  application  of  the  approach 
here  for  Riccati  equations. 

Theorem  4: 

If  Po>0,  A>0,  and  W>0,  then: 

(a)  All  ^>0,  1-0,1,-  •• 

(b)  even,  M odd  with  N<K,  and  M<K. 

(c)  The  following  apriori  bounds  exist  as  N— »CG(or 
as  At — »0)  for  the  continuous  case: 

R0<R2<R,<£R,  C"‘<^lim  PNcC--<iR7<R,<R,<  R, 
46  N-*00or  '331 

<A,t-^0 

with  apriori  bounds: 

Po<po  + A[W+(Por1r1AT</--<C'lim  PN  < •< 

N— *-G3 

<f.  .<f  P +A(W+(P  +A[wr1ATr1r1AT<P  + A[WJ_1AT 
^ u o o 

i.e.  as  many  apriori  bounds  can  be  found  as  is  de- 
sired for  the  numerical  scheme. 

(d)  The  best  estimate  of  PN  (in  the  continuous  case) 
can  be  expressed  as  follows: 

^ - max  (1/2)  [ Rj,  + R(J+1]  (41) 

N 

(e) The  error  in  the  estimate  specified  by  equation 
(41)  can  be  no  worse  than  the  following  bound: 


The  Riccati  Equation 

R 

o 

R, 


Theorem  (4)  contains 


The  proof  of  theorem  (4)  is  outlined  in  Appendix  (C). 


An  iterative  algorithm  is  presented  for  obtaining 
matrix  square  roots.  The  techniques  used  to  derive 
the  algorithm  can  be  used  in  discrete  implementation 
of  the  Riccati  equation.  The  advantage  of  this  ap- 
proach is  the  development  of  apriori  bounds  on  both 
the  matrix  computations  and  the  Riccati  solutions. 
Examples  are  worked  to  illustrate  this  method. 


Appendix  A 


In  this  appendix  all  inequalities  refer  to  matrix 
inequalities,  i.e.  A>B  implies  x^[A-B]x  is  positive 
for  all  vectors  x of  the  appropriate  dimensions.  To 
prove  theorem  2,  lemmas  1-7  must  be  shown  to  be  true 
Lemma  (1): 

Let  A^,  B,  A 2 , , C2  be  positive  definite  matrices 

with: 


A1 

- [ B_1  + C”1 

] -1 

(A.  1) 

a2 

- [ B_1  + C*1 

] ^ 

(A.  2) 

if  then  A^  A . 

Proof : 

A-1  - B-1  + 

1 

cl1 

-1  -1 

-1 

and 

A2  = B + 

C2 

since  Cj^C^ 

and  both  matrices  are 

positive  definite 

-i  -i 

then  we  have 

[17]:C2>  C2 

but 

c-1  - a:1  - 

B-1 

2 2 

-1  -1 

-1 

C1  ' A1  - 

B 

or 

A"1  - B'1  ^ A 

-1 

1 

B_1  (A. 3) 

.-1  . -1 
or  A2  > A| 

, but  since  both 

matrices  are  positive 

definite  implies  A^>  A^. 

q.e.d. 

Lemma  2: 

Let  A,B,  and 

C be  positive  definite 

matrices  with 

C = (A-1  + B' 

'I)”*  , then  the 

following  inequalities 

arise : 


C B (A.  4) 

C CA  (A.  5) 
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Proof  of  Lemma  2: 


-1 


A_1  + B_1 


-1, 


-1 


C >.A 

1^  .,-1 


Note 
Hence 

And  C *J>-B 

Therefore  COi  and  C^A  follows. 

Lemma  3 

The  matrix  subsequences  Rjj  for  N-0,2,4,6,  ‘ ‘ * (all 
even  integers)  form  a monotonically  increasing 
matrix  sequence,  i.e. 


Q.E.D. 


K' 0,<*4<Rfi< 


<%< < 


r2  - I + [if1  +(I+Rq  ) x]  1 

R.  = I + [B_1  +(X+R,)'1]_1 


(A.  6) 
(A. 7) 


But  from  lemma  (1)  if  I+R^^ I+R  , and  using  (A. 6) 
and  (A. 7)  yields: 


R,>  R2 


Now  assume  the  results  hold  for  R^,  it  is  necessary 
to  show  that  they  hold  for  R^+2*  Note  we  assume 

RN^  ^N-2  (N  even) 

but 


Rn  “ 1 + I B_1  +(I+rn_2)"1]"1 


rn+2  - i + [b-1  +(i+RNrY1 

But  from  lemma  (1)  if  (I+R^)2^  I+R^  2 *-mPH-es 
R^+2^  Rjj*  Since  the  results  hold  for  N*0,  and  N*2, 
they  hold  for  all  N.  Q E D 

Lemma  (4) : 

The  matrix  subsequences  R^  for  N*lf3,5,7f"*  (all 

odd  integers)  form  a monotonically  decreasing  matrix 
sequence,  i.e. 


Rx>  R3>  Rj>  Rf>  Rg> 


But  using  lemma  (2)  implies: 


/ S1  - 1 
ni-1  < ~2~ 

. S,+I 


/ X. 

R3  2 


Hence 


R.<Rl 


Sy-I 

For  notation  simplicity  let  B=  > 0 then: 

2 ’ 


-1  ,-l 


Proof : 

The  proof  is  by  mathematical  induction.  By  calcula- 
tion: 


R2  - I+[2(S1-I)_1+(2I)_1]'1>  I - Rq 

v S,-I  . „ 
Hence  R^Rq.  For  notation  simplicity  let  B«_i_  > 0, 
then:  ° ^ 


Now 


R3  - I + [ B +(I+R1)  ] 

R5  “ I + [ B_1  +(I+R3)_1  ]_1 

use  lemma  (1).  If  I+R^^I+R^,  then  this  implies 


from  lemma  (1)  that  Rji*  Rj  or  R^i*  R^  Rj-  Now 
assume  the  results  hold  for  RN>  it  is  necessary  to 
show  that  they  hold  for  Rjj+2*  We  assume: 


But: 


•^<^-2  (N  o^) 


Rn  = i + [ b"1  + (i+Rn.2)'1]'1 


„-l 


1 


RN+2  = I+t  B + U+V  1 

but  from  lemma  (1)  if  I+R^  I+R^_2  this  implies 

RN+2  <V  ®*nce  the  results  hold  for  N+1,3,5,  then 

they  hold  for  all  N.  „ _ „ 

Q.E.D. 

Lemma  (5)  The  following  statements  are  true: 

(A. 5.1)  All  R^  (N  even)  are  bounded  above  by  R^. 

(A. 5. 2)  All  R^  (M  odd)  are  bounded  below  by  Rq. 

(A. 5. 3)  lim  R = R 

N~GO  N N0° 

(A- 5- «)  lim  ^ • R^ 


Proof : 


To  show  (A. 5.1)  is  true,  i.e. 

Ro<R2<R4<R6< <R1  (A-8) 

First  it  is  known  that  all  R^  (N  even)  satisfy: 


N+2 


Sj-i  \-i 


-1,-1 


1 + [l  2 / + (I+V  ] 


Proof 

Again  the  proof  is  by  mathematical  induction.  By 
calculation: 


R3  “ (1/2)  [S1  + I ] 

“(VY1  “| 

' 2 / + I[2I+[2(S1-I)‘1)'1]“1J 


R3  - I + 


srr 


■w  - 1 - ' Vrrv  +(I4Rm)‘1i‘1 

Now  use  lemma  2 which  implies  that: 


^+2 
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or 


^+2 


1 + 1 


(A. 9) 


S,+X 


w,  ■ * 

Since  Rq»I  ^ — j—  for  then  (A.  9)  implies 

that  RN+2<rRl  £or  a^1  N even* 


Q.E.D. 


(for  N odd),  then  from  lemma  (1),  A. -A.  0.  Since 

this  is  true  for  N=3,  it  holds  for1  L all  N. 

Q.E.D. 

Assume  2 ^ or  N even).  By  mathematical 

induction: 

A^  - 0 by  lemma  (1) 

Since  this  is  true  for  N=2,  it  holds  for  all  N. 

Q.E.D. 


To  show  (A. 5. 2)  is  true,  i.e. 

R,i>  R,  "'^R.>  ‘ 


Proof : 


S1  +I.> 


• \ R (M  odd) 

o 

For  N even 

lim 

N-*<X> 

“a  ‘ 

(A.  10) 

For  M odd 

lim 

M— 00 

^ * 

’Vioo 

and 

8 

1 

RM00 

I = R if  S,>  I 
o 1 


"1  2 ' “o  “ “1- 

Hence  R^Jb  R0>  but  all  R^  (M  odd)  satisfy: 


(i+rm)-1  I"1 


‘W  “ 1 + ' (V)  + 

Rjh-2  - 1 - 1 (~r)  + (I+V_1  J"1  > 0 

or  Rm+2-^'  1=r0  f°r  aH  0<ld.  This  implies  (A. 5. 2) 
is  true  for  all  M odd. 

Q.E.D. 

The  statements  (A. 5. 3)  and  (A. 5. 4)  follow  from 
(A. 5.1),  (A. 5. 2),  lemmas  (3,4)  and  the  theorems  of 
monotone  operators  in  a Banach  space  [23]. 

Lemma  (6): 


■W  - **<° 


N+l 


rn>° 


if  N is  odd 


if  N is  even 


Proof : The  following  expression  can  be  calculated: 

■w  - rn  ■ ' (V)  + <i4Vi>'1r1 


-l.-i 


(I+V2>  1 

rewrite  this  as: 

S+l  ' RN  = A1  _ A2 

From  lemma  (1)  we  have:  A^-A^ ^ 0 if  R^  RN-2 

and  A1-A2  <0  if  Rfl_1  RN_2 

but  from  lemma  (5)  we  know  Ro<~R2^Rl 

Hence: 

R,  - R„  > 0 ,N  = 0 

1 o ' 

R2  - Ri  < 0 ,N  = 1 

By  mathematical  Induction  assume  R^_^  <C%_2 


Proof : 

The  proof  is  by  contradiction.  Assume  the  contrary  and 
and 


’Sjod  ” W * 


The  following  expression  is  easily  calculated: 

A = ^JVi  - V “ 

N-^OO  N-^OO 

- [b'1  + a+R^.j)'1]’1 

If  N is  odd,  then  from  lemma  (6)  this  impliesZ^  ^ 0* 
If  N is  even,  then  from  lemma  (6)  this  implies  £>0. 
These  last  two  statements  yield  a contradiction  if 

0.  Note  if7\=  0,  then  everything  is  consistent. 
This  completes  the  proof  of  theorem  (2). 


Appendix  B 


To  show  the  best  estimate  of  a matrix  square  root  is 
is : 

% = max  (1/2)  [ Rjj  + R^  ] 

N 

^ (B.l) 

Define  the  best  estimate  S as  that  estimate  which 
. . . o 

minimizes : 


l/2  1 

min  J = ^ - S ' | (B.  2) 


But  from  lemma  (3)  and  (4)  of  appendix  A we  have 
(for  some  N) : 

Ro^r2<^R4<'  <CR2N<^S1  ^ R2N+1  <v"'<^R3<CR1 

(B.3) 

Now  substituting*?  from  (B.l)  into  (B.2)  yields: 

J - \\l%-  S!1/2)  + 2 <KN-l-Sl/2> 

The  following  inequality  now  results: 

1 , „ .1/2,11  ^ II  1 _ c 1/2 

l(B.4) 


j 46  in  < vsi  mi  + t <vi-sr  > 
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The  inequality  (B.4)  has  the  smallest  upper  bound  if 
and  ^ are  chosen  as  the  largest  values  of  the 

computed  apriori  bounds  in  (B.3).  The  term  best  es- 
timate is  the  guess  of  SQ  which  minimizes  the  right 
hand  side  of  equation  (B.4). 

To  show  that  the  statement  (3b)  of  theorem  (3) 
is  true,  the  worst  case  guesses  of  and  R^i  are 
substituted  into  equation  (B.4).  This  yields  the 
upper  bound  (38b) . 

Q.E.D. 

Appendix  C 

One  would  expect  to  obtain  bounds  on  Riccati  eq- 
uations when  formulated  in  the  factorization  frame- 
work. See,  e.g.  [24*25]  for  simplification  of  Riccati 
equations  when  studied  within  the  context  of  optimal 
control  and  estimation  theory.  The  proof  of  theorem 

(4)  now  follows  as  previously  demonstrated  in  Appen- 
dices A and  B using  the  same  technique.  To  outline 
the  proof  of  theorem  (4)  the  following  lemmas  can  be 
easily  shown  to  be  true  by  following  the  similar 
proofs  in  Appendices  A and  B. 

Lemma  (8) : The  matrix  subsequences  R^  for  N=Q,2,4,6,** 
(all  even  integers)  form  a monotonically  increasing 
matrix  sequence,  i.e. 

Rd<R2^R4<R6<  ••••  < R^  < •••< 

Lemma  (9) : The  matrix  subsequences  R^  for  M=l,3,5,**« 
(all  odd  integers)  form  a monotonically  decreasing 
matrix  sequence,  i.e. 

R^  R3^  R5>  R? 

Lemma  (10) : The  following  statements  are  true: 

(C.10.1)  All  Rjg  (N  even)  are  bounded  above  by  R7  . 
(C.10.2)  All  II.  (M  odd)  are  bounded  below  by  RQ. 

(C.  10.3)  lim  IT  = £ 

(C. 10.4)  lim  % - ^ 


Lemma  (11) 


R.  - E <C0  if  N is  odd 
N+l  N 

R^^  - R^^>  0 if  N is  even 


Lemma  (12)  : For  N even  lim  R^  = R^^ 

N-*00_ 

For  M odd  lim  = 

_ _ M-**CO 

and  =»  R^.  All  statements  in  theorem  (4)  follow 
when  lemmas  (8-12)  are  proved. 
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